Crow-Hassan grassland bird population trends

knitr::opts_chunk$set(echo = TRUE)
library(here)
library(png)
library(jagsUI)
library(ggplot2)
library(ggthemes)
library(tidyverse)
library(kableExtra)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
source(here("Scripts", "5_Prep_Species.R"))

# Set random seed
set.seed(42563)

# Identify species and transsects for analysis.
species_to_run <- prairie_sp
transects_to_sum <- c("CHOA01", "CHOA02", "CHOA09", "CHOA14")

# ID tag for this run and create associated output directory
tag <- "20231128_plus1_to_counts"
mainDir <- here("Results", "991_bird_trends_species")
dir.create(file.path(mainDir, tag), showWarnings = TRUE)
outputDir <- here("Results", "991_bird_trends_species", tag)
subdirectories <- c("annual_n_estimates", "annual_n_plots", "diagnostic_plots", "species_level_estimates", "model_outputs", "synthesized_results", "Rhats", "model_parameters")
lapply(file.path(mainDir, tag, subdirectories), function(x) if (!dir.exists(x)) dir.create(x))

adjust_counts_by_1 <- TRUE

write.table(transects_to_sum, file = here("Results", "991_bird_trends_species", tag, "model_parameters", "transects_summed.txt"))

transect_birds <- bird.tblBirds_Transects

# Merge BB and OA surveys that are actually the same transect
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"


# Filter bird data to only specified transects
transect_birds <- transect_birds %>%
  left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
  filter(SurveyCode %in% transects_to_sum) %>%
  mutate(
    survey_string = paste0(SurveyCode, "-", SurveyDate), # unique code for survey based on area, site #, and date
    year = sub("-.*", "", SurveyDate), # convert survey date to just year
    year = as.numeric(year)
  ) # year to numeric

# Identify days and years in which all transects in the list were actually surveyed
survey_dates <- transect_birds %>%
  group_by(SurveyDate, SurveyCode) %>%
  summarize(survey_count = length(unique(SurveyDate))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(SurveyDate)

survey_complete_years <- transect_birds %>%
  group_by(year, SurveyCode) %>%
  summarize(survey_count = length(unique(year))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(year)

# Re-filter bird data to only dates with complete surveys across all transects. This version forces all transects to have been run on the same day.
transect_birds <- transect_birds %>%
  filter(SurveyDate %in% survey_dates$SurveyDate)

# In years surveyed more than once, pick just the last survey from that year.
surveys_all <- transect_birds %>%
  group_by(year, SurveyCode, survey_string) %>%
  summarize() %>%
  mutate(surveycode_year = paste0(SurveyCode, "-", year))
surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
  ungroup() %>%
  select(survey_string)

# Re-filter bird data based on previous step to only include one survey day each year.
transect_birds <- transect_birds %>%
  filter(survey_string %in% surveys_selected$survey_string)

# Only include counts > 0 (remove rows where bird only detected outside survey area).
transect_birds <- transect_birds %>%
  filter(BirdCountIn > 0)

# Get species list from data and identify species observed at least 5 times.
sp_detected <- unique(transect_birds$BirdCode)
sp_detections <- transect_birds %>%
  group_by(BirdCode) %>%
  filter(BirdCode %in% prairie_sp) %>%
  summarise(
    count_indiv = sum(BirdCountIn),
    n_surveys = n(),
    n_yrs = length(unique(year))
  )
sp_detected_above_threshold <- sp_detections %>%
  filter(count_indiv >= 5) %>%
  pull(BirdCode)

# Re-filter species to run so that it only includes species actually detected in selected surveys
species_to_run <- intersect(species_to_run, sp_detected)

# Also create version with only species detected more than 5 times.
species_to_run_above_threshold <- intersect(species_to_run, sp_detected_above_threshold)

# Save these species list for use in subsequent analyses
saveRDS(sp_detected_above_threshold, here("Results", "991_bird_trends_species", tag, "sp_detected_above_threshold.RDS"))
saveRDS(species_to_run_above_threshold, here("Results", "991_bird_trends_species", tag, "species_to_run_above_threshold.RDS"))

# Function to prepare time series data for species and transects of interest
prep_and_fit_CHBB <- function(species_of_interest = species_to_run, transects_to_assess = transects_to_sum) {
  transect_birds <- bird.tblBirds_Transects

  # Merge BB and OA surveys that are actually the same transect
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"

  # Filter bird data to only those transects
  transect_birds <- transect_birds %>%
    left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
    filter(SurveyCode %in% transects_to_assess) %>%
    mutate(
      survey_string = paste0(SurveyCode, "-", SurveyDate), # Unique code for survey based on area, site #, and date
      year = sub("-.*", "", SurveyDate), # Convert survey date to just year
      year = as.numeric(year)
    )

  # Identify days and years in which all transects in the list were actually surveyed
  survey_dates <- transect_birds %>%
    group_by(SurveyDate, SurveyCode) %>%
    summarize(survey_count = length(unique(SurveyDate))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(SurveyDate)

  survey_complete_years <- transect_birds %>%
    group_by(year, SurveyCode) %>%
    summarize(survey_count = length(unique(year))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(year)

  # Re-filter bird data to only dates with complete surveys across all transects.This version forces all transects to have been run on the same day.
  transect_birds <- transect_birds %>%
    filter(SurveyDate %in% survey_dates$SurveyDate)

  # In years surveyed more than once, pick just the last survey from that year
  surveys_all <- transect_birds %>%
    group_by(year, SurveyCode, survey_string) %>%
    summarize() %>%
    mutate(surveycode_year = paste0(SurveyCode, "-", year))
  surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
    ungroup() %>%
    select(survey_string)

  # Re-filter bird data based on the previous step to only include one survey day each year
  transect_birds <- transect_birds %>%
    filter(survey_string %in% surveys_selected$survey_string)

  # Create df of years with complete data
  survey_years <- unique(transect_birds$year) %>%
    sort() %>%
    as.data.frame()
  survey_years <- rename(survey_years, year = .)

  # Count the number of species of interest across all transects
  sp_counts <- transect_birds %>%
    filter(BirdCode == species_of_interest) %>%
    group_by(year) %>%
    summarize(total_count = sum(BirdCountIn))

  # Zero fill for surveyed years
  sp_counts <- sp_counts %>%
    right_join(survey_years) %>%
    replace_na(list(total_count = 0)) %>%
    arrange(year)

  # This version of state-space model can't be fit when there are zeros in the data (populations that fit zero have no means of recovering/growing). Adding 1 to everything for now...
  sp_counts <- sp_counts %>%
    mutate(
      total_count_orig = total_count
    )

  if (adjust_counts_by_1 == TRUE) {
    sp_counts <- sp_counts %>%
      mutate(total_count = total_count + 1)
  }

  # Identify information about years for model
  yr_min <- min(sp_counts$year)
  yr_max <- max(sp_counts$year)
  year_range <- seq(yr_min, yr_max, 1) # Create a vector of survey years
  n_yrs <- length(year_range) # Counter for how many years of data
  year_range_df <- year_range %>%
    as.data.frame()
  year_range_df <- rename(year_range_df, year = `.`)

  # Add NAs for un-surveyed years
  sp_counts <- sp_counts %>%
    right_join(year_range_df) %>%
    arrange(year)

  #--------------------------

  # JAGS code for Bayesian hierarchical exponential population growth state space model
  {
    sink("Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags")
    cat("
  model {

  # Specify prior values
    init_lnN ~ dnorm(ln1, 1)
    lnN[1] <- init_lnN
    r_0 ~ dnorm(0,1)  # mean intrinsic growth rate
    sigma_r ~ dunif(0, 20) # annual SD in true pop growth rate
    obs_sd ~ dunif(0,20) # annual SD in observation (measurement) error
    # convert SD to precision (tau=1/SD^2)
    sigma_tau <- pow(sigma_r, -2)
    obs_tau <- pow(obs_sd, -2)

  # Likelihood
    for (t in 1:(n_yrs-1)){
      lnN[t+1] ~ dnorm(lnN[t] + r_0, sigma_tau)
    }

  # Observation process
    for (t in 1:n_yrs) {
      lnCount[t] ~ dnorm(lnN[t], obs_tau)
      N_est[t] <- exp(lnN[t]) # convert log(N) back to real scale
    } # close t loop

  # Derived parameters

    # Geometric mean rates of change (%/year between first and last year)
    ch <- N_est[n_yrs] / N_est[1]
    trend <- (ch^(1/(n_yrs-1)))-1

    # Trend from lognormal slope
    slope_trend <- ( n_yrs * sum(1:n_yrs * lnN) - sum(1:n_yrs) * sum(lnN) ) / (n_yrs * sum((1:n_yrs)^2) - sum(1:n_yrs)^2)

  }    # End jags model
", fill = TRUE)
    sink()
  }

  # Bundle data, identify parameters to save, submit to jags
  TRPD_data <- list(
    n_yrs = n_yrs, lnCount = log(sp_counts$total_count),
    ln1 = log(sp_counts$total_count[1])
  )
  parameters <- c("r_0", "sigma_r", "obs_sd", "N_est", "trend", "slope_trend")
  TRPD1 <- jagsUI(TRPD_data,
    inits = NULL, parameters, "Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags",
    n.adapt = 1000, n.chains = 3, n.thin = 10,
    n.iter = 300000, n.burnin = 50000, parallel = TRUE
  )

  #---------------------------------------------

  # Save model output. Models are too big to save everything, not including MCMC chains.
  mod_to_save <- TRPD1[which(sapply(TRPD1, class) != "mcmc.list")]
  saveRDS(mod_to_save, file = paste0(outputDir, "/model_outputs/", species_of_interest, "_model_outputs.RDS"))

  # Assess convergence by verifying all R-hat values < 1.1
  max.rhat <- function(x) {
    max(max_rhat <- c(sapply(x$Rhat, max, na.rm = TRUE)))
  }
  max.rhat.run <- max.rhat(TRPD1) # r.hat <1.1 suggests good convergence
  saveRDS(TRPD1$Rhat, file = paste0(outputDir, "/Rhats/", species_of_interest, "_rhats.RDS"))

  # Summarize r estimates
  r_est <- paste0("r: ", round(as.numeric(TRPD1$mean["r_0"]), 2), " (95% CI: ", round(as.numeric(TRPD1$q2.5["r_0"]), 2), "-", round(as.numeric(TRPD1$q97.5["r_0"]), 2), ")")

  # Summarize percent of r_0 posterior distribution above 0
  f_r0 <- round(as.numeric(TRPD1$f["r_0"]), 2)

  # Data frame with model estimates
  sp_values <- data.frame(
    sp = species_of_interest,
    r0_mean = TRPD1$mean$r_0,
    r0_sd = TRPD1$sd$r_0,
    r0_median = TRPD1$q50$r_0,
    r0_LCI_95 = TRPD1$q2.5$r_0,
    r0_UCI_95 = TRPD1$q97.5$r_0,
    r0_LCI_50 = TRPD1$q25$r_0,
    r0_UCI_50 = TRPD1$q75$r_0,
    f_r0 = round(as.numeric(TRPD1$f["r_0"]), 2),
    sigma_r_mean = TRPD1$mean$sigma_r,
    sigma_r_sd = TRPD1$sd$sigma_r,
    sigma_r_median = TRPD1$q50$sigma_r,
    sigma_r_LCI_95 = TRPD1$q2.5$sigma_r,
    sigma_r_UCI_95 = TRPD1$q97.5$sigma_r,
    sigma_r_LCI_50 = TRPD1$q25$sigma_r,
    sigma_r_UCI_50 = TRPD1$q75$sigma_r,
    f_sigma_r = round(as.numeric(TRPD1$f["sigma_r"]), 2),
    obs_sd_mean = TRPD1$mean$obs_sd,
    obs_sd_sd = TRPD1$sd$obs_sd,
    obs_sd_median = TRPD1$q50$obs_sd,
    obs_sd_LCI_95 = TRPD1$q2.5$obs_sd,
    obs_sd_UCI_95 = TRPD1$q97.5$obs_sd,
    obs_sd_LCI_50 = TRPD1$q25$obs_sd,
    obs_sd_UCI_50 = TRPD1$q75$obs_sd,
    f_obs_sd = round(as.numeric(TRPD1$f["obs_sd"]), 2),
    max.rhat.run = max.rhat(TRPD1),
    method = "SSMr"
  )

  write.csv(sp_values, file = paste0(outputDir, "/species_level_estimates/", species_of_interest, "_species_level_estimates.csv"))

  # Add estimates to simple dataframe
  sp_counts_w_ests <- sp_counts
  sp_counts_w_ests <- sp_counts_w_ests %>%
    mutate(
      N_est_mean = TRPD1$mean$N_est,
      N_est_sd = TRPD1$sd$N_est,
      N_est_median = TRPD1$q50$N_est,
      N_est_LCI_95 = TRPD1$q2.5$N_est,
      N_est_UCI_95 = TRPD1$q97.5$N_est,
      N_est_LCI_50 = TRPD1$q25$N_est,
      N_est_UCI_50 = TRPD1$q75$N_est,
      sp = species_of_interest,
      method = "SSMr"
    )
  write.csv(sp_counts_w_ests, file = paste0(outputDir, "/annual_n_estimates/", species_of_interest, "_annual_N_ests.csv"))

  # Plot posterior draws
  output_path <- here("Results", "991_bird_trends_species", tag, "diagnostic_plots", paste0(species_of_interest, "_diagnostic_plot.png"))
  png(output_path)
  par(mfrow = c(2, 1))

  plot(TRPD1$sims.list$obs_sd, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Observation error", cex = 0.3
  )
  plot(TRPD1$sims.list$r_0, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Process mean", cex = 0.3
  )

  dev.off()

  if (adjust_counts_by_1 == TRUE) {
    sp_counts_w_ests2 <- sp_counts_w_ests %>%
      mutate(
        N_est_LCI_95 = N_est_LCI_95 - 1,
        N_est_UCI_95 = N_est_UCI_95 - 1,
        N_est_mean = N_est_mean - 1,
        total_count = total_count - 1
      ) %>%
      mutate(
        N_est_LCI_95 = ifelse(N_est_LCI_95 < 0, 0, N_est_LCI_95),
        N_est_UCI_95 = ifelse(N_est_UCI_95 < 0, 0, N_est_UCI_95),
        N_est_mean = ifelse(N_est_mean < 0, 0, N_est_mean)
      )
  } else {
    sp_counts_w_ests2 <- sp_counts_w_ests
  }

  # Make time series plot
  plot <- ggplot(sp_counts_w_ests2, aes(x = year)) +
    geom_ribbon(aes(ymin = N_est_LCI_95, ymax = N_est_UCI_95), fill = "grey90", alpha = .6) +
    geom_line(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = total_count), color = "maroon", shape = 1) +
    annotate(geom = "text", label = r_est, x = min(sp_counts_w_ests$year) + 5, y = max(sp_counts_w_ests$total_count), size = 3) +
    scale_x_continuous(breaks = seq(1960, 2021, by = 10)) +
    theme_clean() +
    labs(
      title = paste0(bird.tblBirdSpecies$BirdName[which(bird.tblBirdSpecies$BirdCode == species_of_interest)], " population dynamics"),
      subtitle = paste0("TRPD transects- ", paste(transects_to_assess, sep = "", collapse = ", "), "\nModeled estimate with 95% CI and raw counts \n", r_est),
      x = "Year",
      y = "Total count across transects"
    ) +
    theme(
      plot.title = element_text(size = 11, face = "bold"),
      plot.subtitle = element_text(size = 10, face = "italic"),
      legend.title = element_text(size = 10, face = "plain")
    )

  ggsave(plot = plot, filename = paste0(species_of_interest, "_plot.png"), device = "png", path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), width = 5, height = 3)
}

# # fit state space models for species of interest
# results <- lapply(species_to_run, prep_and_fit_CHBB)
# names(results) <- species_to_run
species_still_to_run <- list.files(
  path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), pattern = NULL, all.files = FALSE,
  full.names = FALSE, recursive = FALSE,
  ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE
)
species_still_to_run <- gsub("\\_plot.png", "", species_still_to_run)
species_still_to_run <- unique(species_still_to_run)
species_still_to_run <- setdiff(species_to_run, species_still_to_run)

library(parallel)
numCores <- detectCores() # get the number of cores available

if (length(species_still_to_run) > 0) {
  lapply(species_still_to_run, prep_and_fit_CHBB)
  # mclapply(species_to_run, prep_and_fit_CHBB, mc.cores = numCores)
}

error_parsing <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  geom_linerange(aes(xmin = sigma_r_LCI_50, xmax = sigma_r_UCI_50)) +
  geom_linerange(aes(ymin = obs_sd_LCI_50, ymax = obs_sd_UCI_50)) +
  geom_label(size = 2.4, label.padding = unit(.075, "lines")) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing

error_parsing_repel <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  ggrepel::geom_label_repel(size = 2.4, box.padding = 0.025, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing_repel

ggsave(plot = error_parsing, filename = "error_parsing_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)
ggsave(plot = error_parsing_repel, filename = "error_parsing_repel_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)

files <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = TRUE)
files_names <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = FALSE)
files_names <- gsub("\\_plot.png", "", files_names)

files <- files[which(files_names %in% species_to_run_above_threshold)]
files_names <- files_names[which(files_names %in% species_to_run_above_threshold)]


for (i in 1:length(files)) {
  cat(paste0("#### ", bird.tblBirdSpecies.ebird$common_name[which(bird.tblBirdSpecies.ebird$BirdCode == files_names[i])], "\n"))
  cat(paste0("![](", files[i], ")\n\n"))
  cat(paste0("**Annual trend (% change/year), 1978-2023:** ", species_level_estimates$r0_mean[which(species_level_estimates$sp == files_names[i])], "% (95% CI: ", species_level_estimates$r0_LCI_95[which(species_level_estimates$sp == files_names[i])], "% to ", species_level_estimates$r0_UCI_95[which(species_level_estimates$sp == files_names[i])], "%)", "\n\n"))
}